/************************************************************************/
/* 
		Auto GGO segmentation

	 Author:				nelfy lee
	 Version:				0.4
	 Modified Data: 2010-04-20
	 Create Date:		2010-04-09
	 E-mail:				nelfylee@gmail.com

	 Changes
	 [0.4/0420]
	  1. improve accuracy.
	 [0.3/0419]
	  1. Fixed Input/Output issues
		2. Removed smoothing output, now all output results are binary
	 [0.2/0413]
		1. add vessel segmentation
*/
/************************************************************************/

#ifndef _NELFY_AUTO_GGO_H_
#define _NELFY_AUTO_GGO_H_


#include "itkObject.h"
#include "itkTimeProbe.h"
#include "itkVectorContainer.h"
#include "itkImageToImageFilter.h"
#include "itkImageRegionConstIterator.h"
#include "itkOtsuThresholdImageFilter.h"
#include "itkConnectedThresholdImageFilter.h"
#include "itkBinaryBallStructuringElement.h"
#include "itkCurvatureFlowImageFilter.h"
#include "itkBinaryMorphologicalClosingImageFilter.h"
#include "itkBinaryMorphologicalOpeningImageFilter.h"
#include "itkBinaryContourImageFilter.h"
#include "itkBinaryDilateImageFilter.h"
#include "itkErodeObjectMorphologyImageFilter.h"
#include "itkObjectMorphologyImageFilter.h"
#include "itkAntiAliasBinaryImageFilter.h"
#include "itkInvertIntensityImageFilter.h"
#include "itkBinaryThresholdImageFilter.h"
#include "itkThresholdSegmentationLevelSetImageFilter.h"
#include "itkFastMarchingImageFilter.h"


namespace itk
{

#define NUMBER_OF_OUTPUTS		10


template<class TInputImage, class TOutputImage>
class ITK_EXPORT AutoGGOSegmentationFilter :
	public ImageToImageFilter<TInputImage, TOutputImage>
{
public:
	typedef AutoGGOSegmentationFilter						Self;
	typedef ImageToImageFilter<TInputImage, TOutputImage>	Superclass;
	typedef SmartPointer<Self>								Pointer;
	typedef SmartPointer<const Self>						ConstPointer;

	itkNewMacro(Self);

	itkTypeMacro(AutoGGOSegmentationFilter,ImageToImageFilter);

	typedef TInputImage                           InputImageType;
	typedef typename InputImageType::Pointer      InputImagePointer;
	typedef typename InputImageType::ConstPointer InputImageConstPointer;
	typedef typename InputImageType::RegionType   InputImageRegionType; 
	typedef typename InputImageType::PixelType    InputImagePixelType; 
	typedef typename InputImageType::IndexType    InputImageIndexType;
	typedef typename InputImageType::SizeType     InputImageSizeType;

	itkStaticConstMacro(ImageDimension, unsigned int, InputImageType::ImageDimension);

	typedef TOutputImage                          OutputImageType;
	typedef typename OutputImageType::Pointer     OutputImagePointer;
	typedef typename OutputImageType::RegionType  OutputImageRegionType; 
	typedef typename OutputImageType::PixelType   OutputImagePixelType; 


	typedef itk::ImageRegionIterator<TInputImage> InputIteratorType;
	typedef itk::ImageRegionIterator<TOutputImage> OutputIteratorType;

	typedef itk::ConstShapedNeighborhoodIterator<TInputImage> InputShapedNeighborhoodConstIteratorType;
	typedef itk::ConstShapedNeighborhoodIterator<TOutputImage> OutputShapedNeighborhoodConstIteratorType;
	typedef itk::ShapedNeighborhoodIterator<TInputImage> InputShapedNeighborhoodIteratorType;
	typedef itk::ShapedNeighborhoodIterator<TOutputImage> OutputShapedNeighborhoodIteratorType;
	typedef itk::NeighborhoodAlgorithm::ImageBoundaryFacesCalculator<TOutputImage> FaceCalculatorType;

	void PrintSelf ( std::ostream& os, Indent indent ) const;

	OutputImageType* GetSmoothedImage()
	{
		return static_cast<OutputImageType*>(this->ProcessObject::GetOutput(1));
	}

	OutputImageType* GetFirstOtsuImage() //2
	{
		return static_cast<OutputImageType*>(this->ProcessObject::GetOutput(2));
	}
	OutputImageType* GetAirRemovedImage()
	{
		return static_cast<OutputImageType*>(this->ProcessObject::GetOutput(3));
	}
	OutputImageType* GetObjectsRemovedImage()
	{
		return static_cast<OutputImageType*>(this->ProcessObject::GetOutput(4));
	}
	OutputImageType* GetLungImage()
	{
		return static_cast<OutputImageType*>(this->ProcessObject::GetOutput(4));
	}
	InputImageType* GetLungWithoutBGImage()
	{
		return lungNoBGImage;
	}
	OutputImageType* GetSecondOtsuImage()
	{
		return static_cast<OutputImageType*>(this->ProcessObject::GetOutput(7));
	}
	InputImageType* GetWeightedImage()
	{
		return weightedImage;
	}
	OutputImageType* GetVesselImage()
	{
		return static_cast<OutputImageType*>(this->ProcessObject::GetOutput(9));
	}
	

protected:
	AutoGGOSegmentationFilter();
	~AutoGGOSegmentationFilter(){}

	InputImagePointer weightedImage;
	InputImagePointer lungNoBGImage;

	void GenerateData();
};




//////////////////////////////////////////////////////////////////////////
// Implementation
//////////////////////////////////////////////////////////////////////////


template <class TInputImage, class TOutputImage>
AutoGGOSegmentationFilter<TInputImage, TOutputImage>
::AutoGGOSegmentationFilter()
{
	this->SetNumberOfRequiredOutputs(NUMBER_OF_OUTPUTS);

	unsigned int numOutputs = static_cast<unsigned int>( this->GetNumberOfOutputs() );
	unsigned int idx;
	if( numOutputs < NUMBER_OF_OUTPUTS )
	{
		for( idx = numOutputs; idx < NUMBER_OF_OUTPUTS; idx++ )
		{
			typename DataObject::Pointer output = this->MakeOutput( idx );
			this->SetNthOutput( idx, output.GetPointer() );
		}
	}
}


template <class TInputImage, class TOutputImage>
void
AutoGGOSegmentationFilter<TInputImage, TOutputImage>
::PrintSelf(std::ostream& os, Indent indent) const
{
	this->Superclass::PrintSelf(os, indent);
	os << indent << "This is auto GGO segmentation" << std::endl;
}



template <class TInputImage, class TOutputImage>
void 
AutoGGOSegmentationFilter<TInputImage,TOutputImage>
::GenerateData()
{
	InputImageConstPointer inputImage = this->GetInput();
	OutputImagePointer outputImage = this->GetOutput();

	OutputImagePointer smoothedImage = this->GetOutput(1);
	OutputImagePointer otsuImage = this->GetOutput(2);
	OutputImagePointer airRemovedImage = this->GetOutput(3);
	OutputImagePointer objectsRemovedImage = this->GetOutput(4);
	InputImagePointer lungImage = InputImageType::New();
	lungNoBGImage = InputImageType::New();
	OutputImagePointer otsu2Image = this->GetOutput(7);
	weightedImage = InputImageType::New();
	OutputImagePointer vesselImage = this->GetOutput(9);

	OutputImagePointer nouse2 = this->GetOutput(5);
	OutputImagePointer nouse3 = this->GetOutput(6);
	OutputImagePointer nouse4 = this->GetOutput(8);



//#define TIMESTART(x)	{ printf(x); itk::TimeProbe clock; clock.Start(); 
//#define TIMESTOP() clock.Stop(); printf("\t\t\t\t\t\t[OK] (time: %.2fs)\n", clock.GetMeanTime());}
#define BG_VALUE					-750		// background value
#define FG_VALUE					100			// foreground value
#define BALL_RADIUS				5				// Rolling ball
#define CANBE_GGO_RATIO		0.85		//


	// image properties
	const unsigned char inputDimension = itkGetStaticConstMacro(ImageDimension);
	unsigned int dim[3] = {1};
	for(unsigned int i=0; i<inputDimension; i++)
	{
		dim[i] = inputImage->GetBufferedRegion().GetSize()[i];
	}

	unsigned long inputImageSize = 1;	
	unsigned long outputImageSize = 1;	
	inputImageSize = dim[0]*dim[1]*dim[2]*sizeof(InputImagePixelType);
	outputImageSize = dim[0]*dim[1]*dim[2]*sizeof(OutputImagePixelType);

	// Create output and internal images

	OutputImageRegionType oregion =  outputImage->GetRequestedRegion();
	InputImageRegionType iregion = inputImage->GetLargestPossibleRegion();

	smoothedImage->SetBufferedRegion(oregion); 
	smoothedImage->Allocate();
	smoothedImage->FillBuffer ( 0);

	outputImage->SetBufferedRegion(oregion); 
	outputImage->Allocate();
	outputImage->FillBuffer ( 0);

	otsuImage->SetBufferedRegion(oregion);
	otsuImage->Allocate();
	otsuImage->FillBuffer( 0);

	airRemovedImage->SetBufferedRegion(oregion);
	airRemovedImage->Allocate();	

	objectsRemovedImage->SetBufferedRegion(oregion);
	objectsRemovedImage->Allocate();	

	lungImage->SetRegions(iregion);
	lungImage->Allocate();	
	lungImage->FillBuffer ( -1024);

	lungNoBGImage->SetRegions(iregion);
	lungNoBGImage->Allocate();	
	lungNoBGImage->FillBuffer ( -1024);

	otsu2Image->SetBufferedRegion(oregion);
	otsu2Image->Allocate();	
	otsu2Image->FillBuffer ( 0);

	weightedImage->SetRegions(iregion);
	weightedImage->Allocate();	
	weightedImage->FillBuffer ( 0);

	vesselImage->SetBufferedRegion(oregion);
	vesselImage->Allocate();	
	vesselImage->FillBuffer ( 0);

	nouse2->SetBufferedRegion(oregion);	nouse2->Allocate();	nouse2->FillBuffer ( 0);
	nouse3->SetBufferedRegion(oregion);	nouse3->Allocate();	nouse3->FillBuffer ( 0);
	nouse4->SetBufferedRegion(oregion);	nouse4->Allocate();	nouse4->FillBuffer ( 0);


	/** 
		1. Smoothing image
		*/

	typedef itk::Image<float,3> ImageFloat3D;
	typedef itk::CurvatureFlowImageFilter<TInputImage, ImageFloat3D> CurvatureFlowImageFilter;

	typename CurvatureFlowImageFilter::Pointer smoothing = CurvatureFlowImageFilter::New();

	//TIMESTART("1. Smoothing...\n");

	smoothing->SetInput(inputImage);
	smoothing->SetTimeStep(0.0625);
	smoothing->SetNumberOfIterations(3);
	smoothing->Update();

	memcpy(smoothedImage->GetBufferPointer(), smoothing->GetOutput()->GetBufferPointer(), outputImageSize);

	//TIMESTOP();
	

	/**
		2. Classify image into two clusters
		*/

	typedef itk::OtsuThresholdImageFilter<ImageFloat3D, TOutputImage>  OtsuThresholdImageFilterType;
	typename OtsuThresholdImageFilterType::Pointer otsu = OtsuThresholdImageFilterType::New();

	//TIMESTART("2. Otsu classifying...\n");

	otsu->SetInput(smoothing->GetOutput());
	otsu->SetOutsideValue(0);
	otsu->SetInsideValue(1);
	otsu->SetNumberOfHistogramBins( 128 );
	otsu->Update();

	memcpy(otsuImage->GetBufferPointer(), otsu->GetOutput()->GetBufferPointer(), outputImageSize);

	//TIMESTOP();


	/**
		3. Remove air
		*/

	typedef itk::ConnectedThresholdImageFilter<TOutputImage, TOutputImage> ConnectedThresholdImageFilter;
	typename ConnectedThresholdImageFilter::Pointer ctFilter = ConnectedThresholdImageFilter::New();

	typedef itk::ImageRegionIterator<TInputImage> InputImageIteratorType;
	typedef itk::ImageRegionIterator<TOutputImage> OutputImageIteratorType;
	OutputImageIteratorType otsuIt(otsu->GetOutput(), otsu->GetOutput()->GetLargestPossibleRegion());
	
	//TIMESTART("3. Removing outside air...\n");

	otsuIt.GoToBegin();
	unsigned int parts = 1;
	unsigned long count = 1;
	for ( unsigned long i = 0; !otsuIt.IsAtEnd(); ++otsuIt, ++i)
	{
		if ( count <= dim[0] || // upper
			count >= dim[0]*(dim[1]-1) || // lower
			count%dim[0] == 1 || count%dim[0] == 0 // border
			)
		{
			if ( otsuIt.Get())
			{
				printf(" [%d] ", parts++);

				ctFilter->SetInput(otsu->GetOutput());
				ctFilter->SetLower( 1);
				ctFilter->SetUpper( 1);
				ctFilter->SetReplaceValue(1);
				ctFilter->SetSeed(otsuIt.GetIndex());
				ctFilter->Update();

				OutputImageIteratorType tmpOrigIt = otsuIt;
				OutputImageIteratorType tmpCTIt(ctFilter->GetOutput(), ctFilter->GetOutput()->GetLargestPossibleRegion());
				tmpOrigIt.GoToBegin();
				tmpCTIt.GoToBegin();

				for ( ; !tmpCTIt.IsAtEnd(); ++tmpCTIt, ++tmpOrigIt)
				{
					if ( tmpCTIt.Get())
						tmpOrigIt.Set(0);
				}
			}
		}

		if ( ++count > dim[0]*dim[1])
			count = 1;
	}

	memcpy(airRemovedImage->GetBufferPointer(), otsu->GetOutput()->GetBufferPointer(), outputImageSize);

	printf("\n");

	//TIMESTOP();


	/**
		4. remove small objects
		*/

	typedef itk::BinaryBallStructuringElement<OutputImagePixelType,3> StructuringElementType;
	typedef itk::BinaryDilateImageFilter<TOutputImage,TOutputImage,StructuringElementType> BinaryDilateImageFilter;
	typedef itk::BinaryErodeImageFilter<TOutputImage,TOutputImage,StructuringElementType> BinaryErodeImageFilter;
	typedef itk::BinaryMorphologicalOpeningImageFilter<TOutputImage,TOutputImage,StructuringElementType> BinaryMorphologicalOpeningImageFilter;
	typedef itk::BinaryMorphologicalClosingImageFilter<TOutputImage,TOutputImage,StructuringElementType> BinaryMorphologicalClosingImageFilter;
	typedef itk::ErodeObjectMorphologyImageFilter<TOutputImage,TOutputImage,StructuringElementType> ErodeObjectMorphologyImageFilter;
	typename BinaryMorphologicalClosingImageFilter::Pointer binMorph = BinaryMorphologicalClosingImageFilter::New();
	typename BinaryDilateImageFilter::Pointer binDilate = BinaryDilateImageFilter::New();

	//TIMESTART("4. Removing small objects...\n");
	
	StructuringElementType structuringElement;
	structuringElement.SetRadius(4);
	structuringElement.CreateStructuringElement();

	binMorph->SetInput(otsu->GetOutput());
	binMorph->SetKernel(structuringElement);
	binMorph->SetForegroundValue(1);
	binMorph->Update();

	memcpy(objectsRemovedImage->GetBufferPointer(), binMorph->GetOutput()->GetBufferPointer(), outputImageSize);

	//TIMESTOP();


	/**
		5. Get Lung
		*/

	typedef itk::ImageRegionConstIterator<TInputImage> ImageConstIteratorType;
	ImageConstIteratorType inputImageIt(inputImage, inputImage->GetLargestPossibleRegion());
	OutputImageIteratorType binImageIt(binMorph->GetOutput(), binMorph->GetOutput()->GetLargestPossibleRegion());
	InputImageIteratorType lungImageIt(lungImage, lungImage->GetLargestPossibleRegion());
	
	//TIMESTART("5. Generating lung image...\n");

	inputImageIt.GoToBegin();
	binImageIt.GoToBegin();
	lungImageIt.GoToBegin();
	for ( ; !inputImageIt.IsAtEnd(); ++inputImageIt, ++binImageIt, ++lungImageIt)
	{
		if ( binImageIt.Get())	// 5.5 Remove background
		{
			lungImageIt.Set(inputImageIt.Get());
		}
	}

	//TIMESTOP();


	/**
		6. Remove background
		*/

	memcpy(lungNoBGImage->GetBufferPointer(), lungImage->GetBufferPointer(), inputImageSize);

	InputImageIteratorType lungNoBGImageIt(lungNoBGImage, lungNoBGImage->GetLargestPossibleRegion());
	
	//TIMESTART("6. Removing lung background...\n");

	for ( lungNoBGImageIt.GoToBegin(); !lungNoBGImageIt.IsAtEnd(); ++lungNoBGImageIt)
	{
		if ( lungNoBGImageIt.Get() < BG_VALUE )
		{
			lungNoBGImageIt.Set(BG_VALUE);
		}

		if ( lungNoBGImageIt.Get() >= FG_VALUE )
		{
			lungNoBGImageIt.Set(BG_VALUE);
		}
	}

	TIMESTOP();


	/**
		7. Get Vessels
		*/


	//TIMESTART("7. Extracting vessels...\n");

	typedef itk::InvertIntensityImageFilter<TOutputImage, TOutputImage> InvertIntensityImageFilter;
	typename InvertIntensityImageFilter::Pointer invert = InvertIntensityImageFilter::New();
	typename BinaryErodeImageFilter::Pointer binErode = BinaryErodeImageFilter::New();
	typename BinaryMorphologicalClosingImageFilter::Pointer closing = BinaryMorphologicalClosingImageFilter::New();

	StructuringElementType structuringElement;
	structuringElement.SetRadius(2);
	structuringElement.CreateStructuringElement();

	binErode->SetInput(objectsRemovedImage);
	binErode->SetKernel(structuringElement);
	binErode->SetErodeValue(1);
	binErode->Update();

	invert->SetInput(airRemovedImage);
	invert->SetMaximum(1);
	invert->Update();

	OutputIteratorType invertIt(invert->GetOutput(), invert->GetOutput()->GetLargestPossibleRegion());
	OutputIteratorType erodeIt(binErode->GetOutput(), binErode->GetOutput()->GetLargestPossibleRegion());
	OutputIteratorType outputIt(vesselImage, vesselImage->GetLargestPossibleRegion());

	invertIt.GoToBegin(); 
	erodeIt.GoToBegin(); 
	outputIt.GoToBegin();

	for (; !outputIt.IsAtEnd(); ++outputIt, ++invertIt, ++erodeIt)
	{
		if ( invertIt.Get() == 1 && erodeIt.Get() == 1)
		{
			outputIt.Set(1);
		}
	}

	structuringElement.SetRadius(2);
	structuringElement.CreateStructuringElement();

	closing->SetInput(vesselImage);
	closing->SetKernel(structuringElement);
	closing->Update();

	memcpy(vesselImage->GetBufferPointer(), closing->GetOutput()->GetBufferPointer(), outputImageSize);

	unsigned int rollingBallRadius = 2;

	typename OutputShapedNeighborhoodConstIteratorType::RadiusType rollingBallRadiusType;
	rollingBallRadiusType.Fill(rollingBallRadius);

	FaceCalculatorType faceCalculator;
	typename FaceCalculatorType::FaceListType faceList;
	typename FaceCalculatorType::FaceListType::iterator fit;

	const float rad = static_cast<float>(rollingBallRadius);

	faceList = faceCalculator( closing->GetOutput(),closing->GetOutput()->GetRequestedRegion(),rollingBallRadiusType);

	for ( fit=faceList.begin(); fit != faceList.end(); ++fit)
	{
		OutputShapedNeighborhoodConstIteratorType closingIt( rollingBallRadiusType, closing->GetOutput(), *fit);
		OutputIteratorType vesselIt(vesselImage, *fit);

		for (float z = -rad; z <= rad; z++)
		{
			for (float y = -rad; y <= rad; y++)
			{
				for (float x = -rad; x <= rad; x++)
				{
					typename OutputShapedNeighborhoodConstIteratorType::OffsetType off;

					float dis = vcl_sqrt( x*x + y*y + z*z );
					if (dis <= rad)
					{
						off[0] = static_cast<int>(x);
						off[1] = static_cast<int>(y);
						off[2] = static_cast<int>(z);
						closingIt.ActivateOffset(off);
					}
				}
			}
		}


		int count = 0;	// how many valid points surround center.
		int index = 0;

		vesselIt.GoToBegin();
		closingIt.GoToBegin();
		for (; !closingIt.IsAtEnd(); ++closingIt, ++vesselIt, ++index)
		{
			if (closingIt.GetCenterPixel())
			{
				typename OutputShapedNeighborhoodConstIteratorType::ConstIterator cit;
				typename OutputShapedNeighborhoodIteratorType::Iterator it;

				count = 0;
				for ( cit = closingIt.Begin(); cit != closingIt.End(); cit++)
				{
					if ( cit.Get() == 1)
					{
						count++;
						if ( count > 6)
							break;
					}
				}

				if ( count <= 5)
				{
					vesselIt.Set(0);
				}
			}
		}
	}

	TIMESTOP();


	/**
		8. Use otsu method again, and then detect the edge of object
		*/

	typedef itk::OtsuThresholdImageFilter<TInputImage, TOutputImage>  OtsuThresholdImageFilterType2;
	typename OtsuThresholdImageFilterType2::Pointer otsu2 = OtsuThresholdImageFilterType2::New();

	//TIMESTART("8. Otsu classifying...\n");

	otsu2->SetInput(lungNoBGImage);
	otsu2->SetOutsideValue(1);
	otsu2->SetInsideValue(0);
	otsu2->SetNumberOfHistogramBins( 64 );
	otsu2->Update();

	memcpy(otsu2Image->GetBufferPointer(), otsu2->GetOutput()->GetBufferPointer(), outputImageSize);

	TIMESTOP();


	/**
		9. Locate GGO positions by comparing otsu2 with outputImage (background removed)
	      1). Construct a ball (radius is 5~8) surround each voxel (value 1) in otsu2.
	      2). Put this ball to outputImage.
	      3). If voxels that inside the ball are mostly vaild (above backgroud value), this voxel is candidate.
		*/


	typedef itk::Vector<int, 4> SeedType;			//value, x, y, z
	typedef itk::VectorContainer<unsigned int, SeedType> SeedContainerType;

	unsigned int seedCount = 0;
	typename SeedContainerType::Pointer seedList = SeedContainerType::New();
	
	//TIMESTART("9. Locating GGO...\n");

	unsigned int rollingBallRadius = 6;

	typename OutputShapedNeighborhoodConstIteratorType::RadiusType rollingBallRadiusType;
	rollingBallRadiusType.Fill(rollingBallRadius);

	FaceCalculatorType faceCalculator;
	typename FaceCalculatorType::FaceListType faceList;
	typename FaceCalculatorType::FaceListType::iterator fit;

	const float rad = static_cast<float>(rollingBallRadius);

	faceList = faceCalculator( otsu2->GetOutput(),otsu2->GetOutput()->GetRequestedRegion(),rollingBallRadiusType);

	InputIteratorType weightedIt;
	

	for ( fit=faceList.begin(); fit != faceList.end(); ++fit)
	{
		OutputShapedNeighborhoodConstIteratorType it( rollingBallRadiusType, otsu2->GetOutput(), *fit);
		InputShapedNeighborhoodConstIteratorType lungnobgit( rollingBallRadiusType, lungNoBGImage, *fit);
		weightedIt = InputIteratorType(weightedImage, *fit);
		
		int shapeCount = 0;

		for (float z = -rad; z <= rad; z++)
		{
			for (float y = -rad; y <= rad; y++)
			{
				for (float x = -rad; x <= rad; x++)
				{
					typename OutputShapedNeighborhoodConstIteratorType::OffsetType off;

					float dis = vcl_sqrt( x*x + y*y + z*z );
					if (dis <= rad)
					{
						off[0] = static_cast<int>(x);
						off[1] = static_cast<int>(y);
						off[2] = static_cast<int>(z);
						it.ActivateOffset(off);
						lungnobgit.ActivateOffset(off);
						shapeCount++;
					}
				}
			}
		}

		
		typename TInputImage::IndexType seedIndex;
		typename OutputShapedNeighborhoodConstIteratorType::OffsetType center;

		center[0] = 0;
		center[1] = 0;
		center[2] = 0;

		int count = 0;	// how many valid points surround center.
		for (it.GoToBegin(), lungnobgit.GoToBegin(); !it.IsAtEnd(); ++it, ++lungnobgit, ++weightedIt)
		{
			if (it.GetCenterPixel())// < 100 && it.GetCenterPixel() > BG_VALUE)
			{
				typename InputShapedNeighborhoodConstIteratorType::ConstIterator ci;

				count = 0;
				for ( ci = lungnobgit.Begin(); ci != lungnobgit.End(); ci++)
				{
					if ( ci.Get() > BG_VALUE)
					{
						count++;
					}
				}

				if ( count > shapeCount * CANBE_GGO_RATIO)
				{
					printf(".");
					weightedIt.Set(count);

					seedIndex = it.GetIndex(center);

					SeedType seed;
					seed[0] = seedIndex[0];
					seed[1] = seedIndex[1];
					seed[2] = seedIndex[2];
					seed[3] = count;

					seedList->InsertElement(seedCount++, seed);
				}
			}
		}
	}

	printf("\n");

	//TIMESTOP();


	/**
		10. Sorting seeds
		*/

	typedef SeedContainerType::Iterator SeedIterator;
	SeedIterator seedIt = seedList->Begin();

	//TIMESTART("10. Processing seeds...\n");

	for ( unsigned int i = 0; i < seedCount; i++, ++seedIt)
	{
		SeedIterator seedIt2 = seedIt;
		seedIt2++;

		for ( unsigned int j = i+1; j < seedCount; j++, ++seedIt2)
		{
			SeedType &iseed = seedIt.Value();
			SeedType &jseed = seedIt2.Value();
			SeedType tmp;

			if ( iseed[3] < jseed[3])
			{
				tmp = iseed;
				
				for ( int m=0; m<4; m++)
					iseed[m] = jseed[m];

				jseed = tmp;
			}
		}
	}

	//TIMESTOP();


	/**
		11. Perform level set segmentation
		*/

	typedef itk::CastImageFilter<TInputImage, ImageFloat3D> CastImageFilter;
	typedef itk::ThresholdSegmentationLevelSetImageFilter< ImageFloat3D, 
		ImageFloat3D> ThresholdSegmentationLevelSetImageFilterType;
	typedef itk::BinaryThresholdImageFilter<ImageFloat3D, TOutputImage>	ThresholdingFilterType;
	typedef  itk::FastMarchingImageFilter<ImageFloat3D,ImageFloat3D> FastMarchingFilterType;


	typename CastImageFilter::Pointer caster = CastImageFilter::New();
	caster->SetInput(lungImage);

	typename ThresholdingFilterType::Pointer thresholder = ThresholdingFilterType::New();

	thresholder->SetLowerThreshold( -1000.0 );
	thresholder->SetUpperThreshold(     0.0 );

	thresholder->SetOutsideValue(0);
	thresholder->SetInsideValue(1);


	FastMarchingFilterType::Pointer  fastMarching = FastMarchingFilterType::New();
	ThresholdSegmentationLevelSetImageFilterType::Pointer thresholdSegmentation =
		ThresholdSegmentationLevelSetImageFilterType::New();

	thresholdSegmentation->SetPropagationScaling( 0.5 );
	thresholdSegmentation->SetCurvatureScaling(15.0);

	thresholdSegmentation->SetMaximumRMSError( 0.02 );
	thresholdSegmentation->SetNumberOfIterations( 2000 );

	thresholdSegmentation->SetLowerThreshold( BG_VALUE);
	thresholdSegmentation->SetUpperThreshold( FG_VALUE);
	thresholdSegmentation->SetIsoSurfaceValue(0.5);

	thresholdSegmentation->SetInput( fastMarching->GetOutput() );
	thresholdSegmentation->SetFeatureImage( caster->GetOutput());
	thresholder->SetInput( thresholdSegmentation->GetOutput() );


	typedef FastMarchingFilterType::NodeContainer           NodeContainer;
	typedef FastMarchingFilterType::NodeType                NodeType;
	NodeContainer::Pointer seeds = NodeContainer::New();

	ImageFloat3D::IndexType  seedPosition;



	//TIMESTART("11. GGO segmenting...\n");

	seedIt = seedList->Begin();
	for ( unsigned int i = 0; i < seedCount; i++, ++seedIt)
	{
		SeedType seed = seedIt.Value();
		printf("\tTrying seed [%d,%d,%d,%d]...\n", seed[0], seed[1], seed[2], seed[3]);


		seedPosition[0] = seed[0];
		seedPosition[1] = seed[1];
		seedPosition[2] = seed[2];

		NodeType node;
		const double seedValue = - 2.0;

		node.SetValue( seedValue );
		node.SetIndex( seedPosition );

		seeds->Initialize();
		seeds->InsertElement( 0, node );

		fastMarching->SetTrialPoints( seeds );
		fastMarching->SetSpeedConstant( 1.0 );

		fastMarching->SetOutputSize(lungImage->GetBufferedRegion().GetSize());

		try
		{
			thresholder->Update();

			// put some conditation here to continue working.
			memcpy(outputImage->GetBufferPointer(), thresholder->GetOutput()->GetBufferPointer(), outputImageSize);
		}
		catch ( itk::ExceptionObject & excep )
		{
			std::cerr << "Exception caught !" << std::endl;
			std::cerr << excep << std::endl;
		}

		break;
	}

	//TIMESTOP();

}


} // namespace itk


#endif


